function forward2()
%
% analyze forward2 data.
% data consists of pulse pairs separated by a variable interval
% analysis includes, as a function of the interval:
% spike counts in both intervals for pulse pairs
% STA analysis of rising phase of second pulse
% including spectral content, slope, variance.

global CONTROL

try
   sf = getmainselection;
   if(sf > 0) 
      pflag = getplotflag;
      QueMessage('Forward2 analysis', 1); % clear the que
      for i = 1:length(sf)
         forward2_2(sf(i),pflag, i);
      end;
   end;
catch
   watchoff;
   QueMessage('Error in Forward2 Analysis routine', 1);
end;


function forward2_2(sf, plot_flag, nth)

global DFILE CONTROL VOLTAGE

persistent h1 h2 h3 h4

err = 0;

dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 1;
% various time definitions:
t_pre = 40; % Time to compute average before a spike; normally 25 msec
t_spec_w = 5; % latest time before spike to measure spectrum (e.g., before foot)
t_post = 15;  % window to display AFTER the spike to show ahp
t_delay = 15; % time after last step before the spike before we consider the data on the rising phase valid

min_isi = 2;

[DFILE, err] = analysis_setup(DFILE, sf);

if(err ~= 0)
   return;
end;

% abstract general information
spike_thresh=number_arg(CONTROL(sf).thresh);
QueMessage(sprintf('Forward2: Spike Threshold: %7.2f mV', spike_thresh));

% information we need to extract
FWD2.dt = []; % the dt array
FWD2.fsl = []; % first spike latency as a function of stimulus
spike_threshf = 0;
time=make_time(DFILE);
tmax=max(max(time));
RATES = (DFILE.rate .* DFILE.nr_channel) / 1000; % array of sampling rates, convert to msec

st = find_steps(DFILE, 'i', 50, 1); % get all step positions in the data set.
recs = length(st);
tx1 = zeros(recs,2);
tx2 = zeros(recs,2);

k = 1;
for i = 1:recs % do data and controls interleaved in this protocol
   if(length(st{i}) == 2) % control has just first steps
      tx1(i,:) = [0 0];
      tx2(i,:) = (st{i})*RATES(i);
   else
      x=st{i}*RATES(i);
      tx1(i,:) = [x(1) x(2)];
      tx2(i,:) = [x(3) x(4)];
   end;
end;
[first_spiket, first_isit, nr_spikest, spike_traint] = find_spikes(DFILE, tx2(:,1), tx2(:,2), spike_threshf);
[first_spike, first_isi, nr_spikes, spike_train] = find_spikes(DFILE, tx1(:,1), tx1(:,2), spike_threshf);

dt = tx2(:,1)-tx1(:,2); % time between pulses (end of first to start of second).
% now get slope before first spike in the second window.
% now calculate slope
for i = 1:recs
   slope(i) = NaN;
   var(i) = NaN;
   if(~isempty(first_spiket(i)) & ~isnan(first_spiket(i)))
      t0 = find_times([tx2(i,1)+10, (first_spiket(i)-2+tx2(i,1))]); % 10 msec delay from onset of step; 2 msec before spike
      if(t0(1) < t0(2))
         t0p = [t0(1):t0(2)]';
         [m, b] = linreg(time(i, t0p), VOLTAGE(i,t0p));
         v2 = VOLTAGE(i,t0p) - (m*time(i,t0p)+b);
         var(i) = sum(v2.^2)/length(v2);
         slope(i) = m;
      end;
      if (max(VOLTAGE(i,t0p))>0)
     %    fprintf(1, 'first: %7.2f tx2: %7.2f Vpk: %7.2f  slope: %7.3f  var: %7.3f, npts: %d\n', first_spiket(i), tx2(i,1), max(VOLTAGE(i,t0p)), slope(i),   var(i), length(v2))
     %    VOLTAGE(i,t0p)
      end;
   end;
end;

fsl = first_spiket;
FWD2.fsl = fsl;
FWD2.dt = dt;
FWD2.slope = slope;
FWD2.var = var;

CONTROL(sf).FWD2 = FWD2;

if plot_flag>=0
   fsize = 7;
   msize = 3;
   msg_pos = [0.37 0.00 0.15 0.07];
   
   if(nth == 1)
      h = findobj('Tag', 'Forward2'); % check for pre-existing window
      if(isempty(h)) % if none, make one
         h = figure('Tag', 'Forward2', ...
            'Units', 'normalized', ...
            'MenuBar', 'none', ...
            'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
            'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
            'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
            'NumberTitle', 'off');
      end
      figure(h); % otherwise, select it
      clf; % always clear the window...
      %Command Menu
      uimenu('Label', 'Close &Window', 'Position', 1, 'callback', 'close(findobj(''Tag'', ''Forward2''));' );
      uimenu('Label', '&Print', 'Callback', 'print;');
      uimenu('Label', 'Print&Setup', 'Callback', 'printdlg;');
      set(h, 'Name', sprintf('Forward2 Analysis - File: %s', CONTROL(sf).filename))
      h1=[]; h2 =[]; h3 = []; h4 = [];
   end;
   
   if(nth > 1)
      h = findobj('Tag', 'Forward2'); % check for pre-existing window
      figure(h);
   end;
   
   if(nth == 1)
      h1=subplot('Position', [0.1 0.5 0.4 0.35]);
   else
      subplot(h1);
      hold on;
   end;
   plot(dt(2:2:end)', first_spiket(1:2:end), '-ko'); % data in black, and thin
   hold on;
   plot(dt(2:2:end)', first_spiket(2:2:end), '-rx'); ...  % overlay curfit in red
      set(gca, 'FontSize', fsize);
   ylabel(CONTROL(sf).V_Unit);
   u=get(gca, 'YLim');
   %set(gca,'XTick',[ ]);
   v=get(gca, 'XLim');
   v(2)=tmax;
   set(gca, 'Xlim', v);
   datac_setcrosshair(gca, 'Forward2_info', 'dt(ms)', 'fsl(ms)', msg_pos);
   
   if(nth == 1)
      h2 = subplot('Position', [0.1 0.1 0.4 0.35]);
   else
      subplot(h2);
      hold on;
   end;
   plot(dt(2:2:end), first_spiket(2:2:end)-first_spiket(1:2:end), '-ko'); % data in black, and thin
   set(gca, 'FontSize', fsize);
   ylabel(CONTROL(sf).V_Unit);
   u=get(gca, 'YLim');
   %set(gca,'XTick',[ ]);
   v=get(gca, 'XLim');
   v(2)=tmax;
   set(gca, 'Xlim', v);
   datac_setcrosshair(gca, 'Forward2_info', 'dt(ms)', 'deltafsl(ms)', msg_pos);
   
   h3 = subplot('Position', [0.55 0.5 0.4 0.18]);
   [cs, vcs] = mean_var(slope(1:2:end));
   vcs = sqrt(vcs);
   plot(dt([2,end]), [cs cs], 'k-');
   hold on;
   plot(dt([2,end]), [cs+vcs cs+vcs], 'k--');
   plot(dt([2,end]), [cs-vcs cs-vcs], 'k--');
   plot(dt(2:2:end), slope(1:2:end), 'ko', 'markersize', 2);
   plot(dt(2:2:end), slope(2:2:end), 'rx');   
   datac_setcrosshair(gca, 'Forward2_info', 'dt(ms)', 'slope(mV/ms)', msg_pos);
   
   h3 = subplot('Position', [0.55 0.70 0.4 0.18]);
   [vs, vvs] = mean_var(var(1:2:end));
   vvs = sqrt(vvs);
   plot(dt([2,end]), [vs vs], 'k-');
   hold on;
   plot(dt([2,end]), [vs+vvs vs+vvs], 'k--');
   plot(dt([2,end]), [vs-vvs vs-vvs], 'k--');
   plot(dt(2:2:end), var(1:2:end), 'ko', 'markersize', 2);
   plot(dt(2:2:end), var(2:2:end), 'rx');   
   datac_setcrosshair(gca, 'Forward2_info', 'dt(ms)', 'var(mV^2/ms)', msg_pos);
   
   h4 = subplot('Position', [0.55, 0.1, 0.4, 0.35]);
   tx = find_times([tx2(1,1) tx2(1,2)]);
   for i = 1:recs
      tx = find_times([tx2(i,1) tx2(i,2) (first_spiket(i)-2+tx2(i,1))]);
      if(isnan(tx(3)))
         tx(3) = tx(2);
      end;
      txp{i}=[tx(1):tx(3)];
      
      plot(time(i,txp{i})', VOLTAGE(i,txp{i})');
      hold on;
   end;
   datac_setcrosshair(gca, 'Forward2_info', 'T(ms)', 'V(mV)', msg_pos);
   
end;
